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Taylor-Couette flow in the presence of a magnetic field is a problem belonging to classical 
hydromagnetics and deserves to be more widely studied than it has been to date. In the 
nonlinear regime the literature is scarce. We develop a formulation suitable for solution 
of the full three dimensional nonlinear hydromagnetic equations in cylindrical geometry, 
which is motived by the formulation for the magnetic field. It is suitable for study at finite 
Prandtl numbers and in the small Prandtl number limit, relevant to laboratory liquid 
metals. The method is used to determine the onset of axisymmetric Taylor vortices, and 
finite amplitude solutions. Our results compare well with existing linear and nonlinear 
hydrodynamic calculations and with hydromagnetic experiments. 



1. Introduction 

The motion of an incompressible viscous fluid between concentric rotating cylinders is 
one of the most important problems of fluid dynamics and i s much studied as a benchmark 
to investigate issue of instability and nonlinear behaviour. Taylor (1923)| found that if the 
rotation of the inner cylinder is greater than some critical value then circular-Couette 
flow (CCF) becomes unstable to axisymmetric perturbations. A secondary flow appears 
which has axial and radial motion in the form of pairs of toroidal vortices, now known 
as the Taylor-vortex flow (TVF). If the inner cylinder is driven further then this flow 
becomes unstable to non-axisymmetric perturbations. Azimuthal waves appear in the 
Taylor-vortices and the whole pattern rotates at some wavespeed (wavy modes). 

In his landmark 1961 book on stability theory, Chandrasekhar devoted equal attention 
to the hydrodynamic and the hydromagnetic Couette problems; the latter is the case in 
which the fluid is a conducting liquid (e.g. mercury, liquid gallium, liquid sodium) and 
a magnetic field is applied externally. Despite this early interest in the hydromagnetic 
Couette problem, which included experiments performed by Donnelly & Ozima (1962) 
and by Donnelly & Caldwell (1963), most of the activity of the following years was devoted 
to the hydrodynamic case. Among the few studies of the effects of the magnetic field it 
is worth remembering the works by Velikhov (1959)| , [Kurzweg (1963) , Roberts (1964) 



who extended Chandrasekhar's theory to non-axisymmetric bifurcations from circular- 



Couette flow, Chang & Sartory (1967) and Baylis & Hunt (1971) at finite aspect ratio. 



Later Tabeling (1981), using a method similar to Davey's (1962) amplitude expansion, 
calculated effective viscosity of axisymmetric flow in the Taylor vortex flow regime; he 
compared against Donnelly's (1962) experiments which indicate that the onset of wavy 



vortices is significantly inhibited by the magnetic field. Nagata (1996)| has more recently 
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investigated nonlinear solutions in the planar geometry, and Hollcrbach (2000a) shows 
Taylor cells in spherical geometry 

The aim of this paper is to investigate effects induced on Couette flow by an externally 
applied magnetic field. This paper is meant to be the first of a series and is dedicated to 
the development of a suitable formulation for solving numerically the governing nonlin- 
ear three dimensional magnetohydrodynamic (MHD) equations in the cylindrical Couette 
geometry. The numerical method which we propose can be used for any value of radius ra- 
tio, is suitable for time stepping, has good stability features, is relatively easy to program 
and is more accurate than existing methods. 

Our work is also motivated by the renewed interest in MHD flows in confined geome- 
tries which arises from current and planned experiments to produce dynamo action in 
the laboratory (Gailitis et al. 2001; Stieglitz & Muller 2001). It must be stressed that our 
work does not apply directly to the dynamo problem for two reasons. Firstly the cylin- 
drical Couette configuration is a possible geometry for these studies, but it has not been 
used in experiments yet, and it may prove not to be not the most efficient one (Laurc, 
Chossat & Daviaud, 2000). Secondly, our investigation refers primarily to bifurcations at 
relatively small Reynolds numbers, while the dynamo experiments require rather large 
Reynolds numbers due to the small magnetic Prandtl number of liquid metals. Despite 
these two limitations, our work is related to the dynamo problem because it is important 
to have precise results at small Reynolds number MHD flows in order to develop and test 
modern acoustic flow visualisation techniques, Kikura, Takeda & Durst (1999), which 
offer the best chance to detect flow patterns in MHD dynamos. The lack of flow visual- 
isation has clearly held back progress in the hydromagnetic Couette problem compared 
to the hydrodynamic case. 



2. Equations 

The equations governing incompressible hydromagnetic flow are 

d t u+(u- V)u = --Vp + vV 2 u+ — CVAB) AB, V • u = 0, (2.1a, b) 
P PPa 

d t B = \V 2 B + V A (it A B). V • B = 0, (2.1c, d) 

where u is the fluid's velocity, B the magnetic field, p the pressure, p the density, v 
the kinematic viscosity, A the magnetic diffusivity and the magnetic permeability. 
Hereafter we assume that p, v, X and po are constant. The fluid is contained between 
two concentric cylinders of inner radius i?i and outer radius i?2- The inner and outer 
cylinder rotate at constant angular velocities f2i and fl 2 respectively. A magnetic field 
Bq = p.oHz is applied externally in the axial direction. We make the usual simplifying 
assumption that the cylinders have infinite length and use cylindrical coordinates (r, 9, z). 

Throughout the rest of this work we will make the variables dimensionless using the 
following scales: 

5 = i?2 — Ri, length (gap width); 8 2 /v, time (viscous diffusion); 
v/5, velocity; poH magnetic field. 

We introduce the following dimensionless parameters: radius ratio (tj), Reynolds numbers 
(Re\ and Re2), Hartmann number (Q) and magnetic Prandtl number (£) defined as 



Ri^iS . _ n _ p 2 H 2 aS 2 



v 



r, = R 1 /R 2 ; Re, = ^^, t = 1, 2; Q = ^ ; £ = -. (2.2) 

v pv A 
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The dimensionless forms of (2.1 a, c) are then 



d t u + (u ■ V)u = -Vp + V 2 m AB) AB, (2.3a) 

dtB = )^ 2 B + V A (it A B). (2.3b) 

A steady-state solution of the governing equations is circular-Couette flow, u — (0, u$, 0) 
where ug is ug = Ar + B/r. The constants A and B are determined by the no-slip 
boundary conditions. We set u = u + u' and p = p + p' . The deviation, u' , then 
satisfies the homogeneous Dirichlet bou nda ry condition, u' = at R\ , i?2 . Subtracting 
the Navier-Stokes equation for u from (2^ a), the evolution of u' is now described by 

(d t - V 2 )u' = N - Vp', V • u' = 0, (2.4a, 6) 

(0 t - iv 2 )S = 7V fl , V- J B = 0, (2.4c, d) 

with nonlinear terms, 

JV=|(VAB)AB-(«-VX-(«'-V)ii, N B = V A (it A B). (2.4e,/) 

The magnetic Prandtl number £ is very small in liquid metals available in the lab- 
oratory, so we set B = B + £,b. In the limit £ — > the governing equations become 

(ft - V 2 )tt' = AT - Vp', V • «' = 0, (2.5a, 6) 

V 2 b = N B , V-6 = 0, (2.5c, d) 

where, 

N = Q(V A 6) A Bo - {u ■ V)u' - (u' ■ V)u, N B = -V A (u A B a ). (2.5e,/) 
Note that these equations are descriptive rather than predictive for b. 



3. Boundary conditions 



The governing equations (2.4) represent a tenth order system in r and we therefore 
require ten boundary conditions. The first six are simply the no-slip condition, u r — 
ug — u z — applied at the boundaries r = R\ and r = R%. The boundary conditions 



for the magnetic field depend on the nature of the container, as discussed by Roberts 



(1964), who determined conditions for arbitrary values of electrical conductivity. The 



experiments by Donnelly & Ozima (1962) used mercury with Perspex and stainless-steel 
containers. Only a small difference was found between the results obtained using different 
containers. Hereafter we consider only the simple case of insulating boundaries. 

Ampere's law says that, J = £ -1 V A B ~ Q, when r < i? x or r > R 2 , as the current 
within an insulator must be zero. It follows that the magnetic field is irrotational and 
can be expressed in terms of a potential, tjj, in the following way: 

B = -Vi/j, -V • B = V 2 V> = 0. (3.1) 

This equation can be solved for tp by separation of variables, ip(r, 8, z) = R(r) 0(9) Z(z). 
In our periodic coordinates we obtain 



(3.2) 
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where m is integer. The equation for R(r) satisfies the modified Bessel equation, 



i R'(r) + R"(r) - (a 2 + ^ R(r) = 0. 



(3.3) 



The boundary conditions for R(r) then depend on the type of solution. 

If ip is independent of 9 and z (m — a — 0) then tp must be constant and so B = 0. 
But this means that we have three conditions at each boundary, and we only need two. 
However, the divergence-free condition implies that a solution which is independent of 9 
and z must have no radial component. It is therefore sufficient to take 

B e = B z = 0. (3.4) 

If tp is independent of z but depends on 9 (a = 0, m ^ 0) then R(r) — r . Recalling 
that B = V-0, we have 

d e B r = ±mB e , B z = 0. (3.5) 

If %p is z dependent (a ^ 0) then R(r) = B m (r) where B m {r) denotes either of the 
modified Bessel functions I m (ar), K m (ar). We obtain 

d z B r = f& B z , - d e B z = d z B e . (3.6) 
B m r 

In the outer region r > R 2 the field tends to zero, B — > as r — ► oo, and in the inner 
region r < R\ B must remain finite, which implies that we take 

m s / I m (ar) r+ m r < R l 

R(r) = < i \ or _ m on (3.7) 
v ' I K m (ar) r m r > R 2 



The two relations (3.4), (B.5), or (3.6 ) applied at the points r = R\ and r = R 2 , given 
the appropriate function from (3.7), are equivalent to Roberts' insulating boundary con- 
ditions. In this way we have the remaining four boundary conditions which are required. 



4. Formulation and solution 

The difficulty with primitive-variable formulations is how to ensure a divergence-free 
field. A possible solution is to combine time splitting and pressure projection. The diver- 
gence of the momentum equation gives a Poisson equation for the pressure which is used 
to project the velocity into the space of solenoidal functions. No pressure term occurs 
naturally in the induction equation, and if not removed divergence can build up in the 
solution for the magnetic field, especially at larger magnetic Prandtl numbers. However, 
an arbitrary projection function could be added in order for the divergence to scale with 
the timestep. Marcus (1984) used an influence matrix method in order to implement the 
correct boundary conditions for the pressure (Rempfer, 2002). This technique leads to 
a divergence that is zero to machine accuracy. In both methods an adjustment to the 
field is made at each timestep. However, for the hydromagnetic case there is no timestep 
in the small Prandtl number limit; without potentials it is difficult to invert the Pois- 
son equation for the magnetic field whilst simultaneously ensuring it is divergence- free. 
We propose a formulation able to cope with both finite Prandtl numbers and the small 
Prandtl number limit without significant adjustments. 

Popular in MHD is the toroidal-poloidal potentials form where variables are decom- 
posed as A = V A (ipe) + V A V A (4>e) where e is a vector constant. To eliminate the 
pressure in the momentum equation it is commonplace to take the e-components of the 
first and second curls as the governing equations for the velocity. For the magnetic field 
it is sufficient to take the e-components of the induction equation and its first curl. 
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In spherical geometry one assumes e = r s , the spherical radius. Taking r s -components 
of successive curls, the Naiver-Stokes and induction equations separate into one equation 
for each of the potentials, Hollerbach (2000&). Unfortunately complications can arise if 
this method is used in cylindrical geometry. The choice e = z leads to separate equations 



for each of the potentials but raises the order of the equations in r. However, Marques 



(1990) has derived the extra boundary conditions required for the hydrodynamic problem. 

Taking the second curl leads to an operator acting on one of the potentials in the form 
of double Laplacians. In spherical geometry Tilgner & Busse (1997) successfully imple- 
mented a second order code using this formulation with stress- free boundaries. Hollerbach 
(2000&) also used this formulation with no-slip boundary conditions but found it to be 
unstable, even for very small timesteps, unless an implicit first order time discretisation 



was used. In the cylindrical geometry Rudiger & Feudel (2000) used the formulation of 
Marques (1990), but similarly used a first order method to avoid numerical difficulties. 
It is not at all clear that higher derivatives necessarily entail numerical instability in 
hydrodynamical solvers. However, for reasons that become apparent when the velocity 



field is discussed (§4.3), the second curl is avoided. 

Instead we are motivated by a parallelism with the magnetic field and take only the 
first curl. As the pressure has not been eliminated, we also take the divergence of the 
momentum equation. A single-curl formulation was proposed by Glatzmaier (1984) in 
the spherical geometry and has proved very successful. 

We have made the choice e = r where r is the cylindrical-polar radius. This choice gives 
equations which couple the potentials, but this is no particular problem; with primitive 
variables r and 9 components are coupled. Fortunately, although we still raise the order 
of the equations, the extra derivatives appear in the periodic coordinates and so no extra 
boundary conditions are required. To ensure capture of all possible solutions extra terms 
along 8 and z are added, in order to accommodate solutions that are independent of 
both 9 and z. Full expansion of variables has the form 



A = t/jQ 8 + 0o z + V A (tpr) + V A V A (<jyr), 



(4.1) 



where ip(r,t,z), <j)(r,t,z) and ?/>o( r )j 4>o{ r ) contain the periodic and non-periodic parts 
respectively. We discuss first the formulation for the magnetic field, as it motivates the 
method for the velocity. 



4.1. The magnetic field 

The magnetic field is expanded as 

B = T d + V z + V A{Tr) + V AV A(Vr), (4.2) 

and substituted into the induction equation, ( |2.4| rf). For the non-periodic potentials To, Vo 
the governing equations are obtained from the 6, z components 



{d t -\{V 2 --))% = §■ Ni 
? r 

(d t ~^W 2 )P = z-N B . 
with boundary conditions at R2, 

T = 0, V = 0. 



(4.3a) 
(4.36) 

(4.4) 



The periodic potentials T,V are assumed to be of form e 1 ( Q2; + me ) j n order to match the 
boundary conditions a spectral expansion will be required. There is no pressure term to 
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eliminate here, so we take the r-components of the induction equation and its first curl 

9 1 - _ 9 1 

P = — T--iV B , (4.5a 



-W 2 c (d t - - V 2 )T - — d ree T +^(d t - -V 2 )d ez V ^-r-VANi 



where 



£r 3 



v 2 



+ - ft., 

r 



(4.56) 
(4.6) 



with boundary conditions at R\ , i?2 
a = : 



a e r = o, ( v* ± — ft. ) ft>p 



a ^ 



v 2 r- 

-d e T 

r 



■ d 0z V = 0, 



B„ 



8 r Br 



- + ft. 5^ 
r 



0. 



(4.7a) 
(4.76) 



This formulation is suitable for the small Prandtl number limit (2.5), relevant to lab- 
oratory liquid metals. We expand b by the same potentials and make the replacements, 



ft ->0, 



1. 



throughout. The field b satisfies the same boundary conditions as B. 

Having settled on a formulation of the equations, we now discuss the numerical method. 

4.2. Numerical method 

It is customary to adjust the radial range into the unit interval, so a new radial variable 
x is defined as 

r = Ri + x, Ri = 77/(1-77), ire [0,1]. (4.8) 

If a field is expected to have mi-fold rotational symmetry, such as the case of wavy 
modes, then variables are expanded as 



JY 



A{x,6,z,t) = Y J E E A nkm (t)T*(x)e^ kz+m - me 

n=0 |fc| <K |m| <M 



(4.9) 



on the domain [0,1] x [0, 2ir/m{\ x [0,27r/a] where T*(x) is the n th shifted Chebyshev 
polynomial. Variables are collocated on the N + 1 extrema of Tjy(x). This arrangement 
of points is well suited to our problem with the points concentrated near the boundaries. 

As the velocity and magnetic field are coupled by the nonlinear terms it makes sense to 
treat them explicitly. Nonlinear terms are evaluated pseudospectrally, where necessary. 
Large terms in the zero-modes, like circular-Couette flow and the imposed magnetic field, 
can be extracted and calculated exactly. Let q indicate the time discretisation t q — q At 
with q = 0, 1, 2, .... We choose to use second order Adams-Bashforth to estimate Nb at 
the intermediate time q + h . 

The linear terms are easier to evaluate and are timestepped using the implicit Crank- 
Nicolson method. Substituting the spectral expansion in the governing equations (4.5) 
and the boundary conditions (4.7), after collocation the problem for each k,m mode 
becomes 



X 



T 
V 



9+1 



= Y 



' T ' 


1 




V 


+ 


. ^2 _ 



(4.10) 
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where X and Yare matrices. The vector [T, V\ q contains the spectral coefficients of the 
potentials at time t q . The nonlinear terms have been evaluated on the collocation points. 
Despite the Fourier expansions, the matrices X, Y are real and are calculated by the same 
routine, as they differ only by a scalar constant, namely ±At. 



With this formulation modifications for the small Prandtl number limit (2.5) are rel- 
atively minor. The magnetic field is now completely defined by the velocity at some 
particular time t q . For each k, m mode the problem in matrix- vector form is now simply 

(4.11) 



X 



' T ' 


q 




V 







The matrices X are calculated by the same routines above, as, again they only differ by 
scalar constants. 

4.3. The velocity field 

Unless there is an externally imposed velocity field, there is no non-periodic pressure to 
be eliminated. The non-periodic part of the velocity is then treated in the same manner 
as the magnetic field, i.e. 

(d t ~ (V 2 — -))ipo = 8 • N, (4.12a) 



(d t - V 2 )0 O = z ■ N, (4.126) 

with boundary conditions 

^o = 0o = O, (4.13) 

on Ri , i?2 ■ 

For the periodic part we follow the procedure applied to the magnetic field and take 
only the first curl. As the pressure has not been eliminated, we also take the divergence 
of the momentum equation. The equations obtained are 



~2 d0z ^ 
-V 2 (d t -V 2 )V> 



V 2 (5 t 
2 



2 1 

— d r8o4>= r ■ (N - Vp), 



r* 

72., 



:{d t 



2V 2 )d e 
N. 



1 



r ■ V A AT, 



(4.14a) 

(4.146) 
(4.14c) 



These equations may look compicated enough, but taking the second curl they are much 
worse! However, they simplify considerably for the axisymmetric problem. N ote also the 
s impl ification that the linear differential operators on the left-hand sides of (4.5 a, 6) and 
(4. 14 a, 6) are the same with £ — > 1. Fourth order derivatives have been avoided, otherwise 
as |d p T*(x)/dx p \ = 0(n 2p ) matrices can become difficult to invert accurately with larger 
truncations. 

Every governing equation is only second order in r, and therefore all equations have 
the same number of associated boundary conditions. This permits us to take the same 
radial truncation N for all variables, so all matrices are likewise of the same size. This 
simplifies the actual implementation enormously! In fact, we timestep the governing 



equations (4. 14 a, 6), the same way as the magnetic field, 



X 



- i, - 


= y 




9 

+ 






. $ . 




. ^2 . 



9+ 



(4.15) 



and the matrices X , X 1 Ymay be precomputed by the routine for the magnetic field as 
the equations are the same, but for scalar constants. Thus, linear terms for the velocity 
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are also timcsteppcd using Crank-Nicolson. Using this method, Marcus (1984) found that 
a numerical neutrally stable oscillation can occur at large wavenumbers. Fortunately this 
does not present a problem here, as the presence of the magnetic field tends to reduce 
the natural wavenumber. 

Together, the evolution equations (i.l4a,b) are only fourth order in r for the poten- 
tials, but there are six boundary conditions for the velocity. They are timestepped with 
boundary conditions ug = u z = 0, or 



rd z ip + d re <j> = 0, -dei> + (2 + rd r )d z <p = 0. 



(4.16a, 6) 



The pressure-Poisson equation (4.14 c) is inverted with the boundary condition u r = 0, 
or, 



0, 



(4.16c) 

essentially the no-penetration condition. 

Inversion of ( 4.1 4| c) requires a boundary condition, indirectly determined by ( 1.16j c), 
in terms of p. The adjustment for pressure is 

9+1 



x- 







(4.17) 



In order to separate (f> q+1 from ip q+1 we must work with X 1 rather than X. For each 
k, m mode "S7 2 is just a scalar. Imposing —rV 2 4> q+1 = the boundary condition becomes, 



X 



-d r p) = 



(4.18) 



where X is the lower left quadrant of X 1 . Condition ( 4.1Sj ) is implemented in the 
usual manner by multiplying on the left with the coefficients T*(x) at the boundaries. 



5. Results 

5.1. Linear stability 

The linear part of the code is shared between the velocity and magnetic fields. Appropri- 
ate tests are determining the critical Reynolds number, Re c , for the onset of Taylor-vortex 
flow in the presence/absence of a magnetic field, and determining the growth/decay rates 
of cither field. 

An eigenfunction of the linearised equations grows or decays exponentially at a rate a. 



Barenghi (1991) examined convergence with At for the velocity by comparing against a 
known growth rate. A simple initial disturbance to the appropriate mode is <fi oc x 2 (l — 
x) 2 smaz, or equivalently 4>o,±i,o = 3A, <p2.±i,o = — 4A, </>4,±i,o = A, which satisfies the 
boundary conditions and mimics TVF surprisingly well. 

To ensure that the boundary conditions for the magnetic field have been set up correctly 
we check our method against analytically derived decay rates (see the appendix). Table 
|l| shows results of the test of growth rates and the comparison with Barenghi (1991). 
Note that the error is proportional to At 2 . 

To check the interaction of the two fields and the case £ — * we compare the onset of 
TVF against Roberts (1964). Table || shows the number of modes required to reproduce 



a few of Roberts' results to five significant figures. 

5.2. Nonlinear two-dimensional flow 

The saturation to a steady flow, for some not too large Re\ > Re c , provides a testing 
ground for the evaluation of nonlinear terms. Barenghi (1991)| compares values for veloc- 
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At 


% error in 


a 


% error in as 


0.01 


5.74x10" 


12 




0.003 


5.18x10" 


3 


1.48xl0~ 2 


0.001 


5.76x10" 


4 


1.65xl0 -3 


0.0003 


5.2x10" 


5 


1.5xl0" 4 


0.0001 


6x10" 


6 


2xl0" 5 



Table 1. Error in growth and decay rates. N — > oo. For rj = 1/1.444, a = 3.13, 7?ei = 80, 
Re2 — the growth rate is a = 0.430108693 (Barenghi 1991). For the magnetic field m = 1, 
<7s£ = 14.055585, f = 1. The error is proportional to At 2 . 



Q 




N 


Re c 


30 


2.69 


8 


280.97 






10 


281.05 


100 


1.73 


8 


463.20 






10 


463.52 


300 


0.928 


8 


796.52 






10 


798.57 






12 


798.55 



Table 2. Critical Reynolds numbers for varying numbers of modes and magnetic field strengths. 
rj = 0.95 with insulating walls. For the largest number of modes in each case the values are the 
same to five significant figures as the results of Roberts' (1964) calculations. 



N 


K 








10 


6 


4.236577 


17.94932 


33.48869 




8 


4.236615 


17.97669 


33.66495 




12 


4.236616 


17.97902 


33.70222 


16 


6 


4.233596 


17.94086 


33.45839 




8 


4.233635 


17.96816 


33.64135 




12 


4.233635 


17.97046 


33.67982 






4.23363 


17.9705 


33.6805 



Table 3. Radial velocity at the outflow. At — » 0. rj = 0.5, x — 0.5, a — 3.1631, and a fixed 
outer cylinder. In the last row are the values for Rei = 72.4569, 106.066, 150.000 (Jones 1985). 



ities at the outflow which are in agreement with results obtained by Jones (1985) using a 
different method. In table |^ we examine convergence with N,K. Generally we find that 
convergence is quicker in z than r, and, for a given truncation accuracy decreases with 
increasing Re\. More energy is found in the higher modes as Re± is increased. 

5.3. Wavespeeds in wavy TVF 

A simple small wavy perturbation to axisymmetric TVF that satifies the boundary con- 
ditions is ip oc x 2 (l — x) 2 sinmi(9 for some wavy mode m\. The purtubation will either 
decay or grow and saturate depending on whether or not the parameters are in the wavy 



TVF regime. King et al. (1984) compared wavespeeds founds from physical experiments 
with numerical calculations. They found that "the wavespeed is a sensitive indicator of 
the accuracy of a numerical code" . . . "any compromise in numerical resolution changes 
the wavespeed by several percent" . They also argue that the wavespeed can be measured 
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Rei/Re c 2n/a Marcus Measured 

3.98 2.40 0.3443 ± 0.0001 0.3440 ± 0.0008 
3.98 3.00 0.3344 ± 0.0001 0.3347 ± 0.0007 
5.97 2.20 0.3370 ± 0.0001 0.3370 ± 0.0002 

Table 4. Wavespeeds expressed as a fraction of the angular velocity of the inner cyliner. 

rj = 0.868, Re c = 115.1, mi = 6. 



more precisely, both in experiment and numerically, than torques which are dependent 
on axial wavelength (see the c omparison with torque experiments in §|J). 

A few of the results used by Marcus (1984) as a test for his numerical method are given 
in table ^. Marcus' numerical results were well within the range of experimental error of 
about 1%. Calculations with our code gave results all within 0.1% of Marcus' values. 



6. Nonlinear hydromagnetic flow and comparison with experiment 

It should first be noted that the results of the previous section agree well with ex- 
periments. In this section we directly compare our results with hydromagnetic torque 
experiments. 

The torque per unit axial length on the inner cylinder is defined as 

,2tt ,27r/a / x \ 

G=— rd6 dz(--d r )ug . (6.1) 

2 W Jo \r J 

r=R 1 

There is no magnetic torque with insulating boundaries. For an axisymmetric flow and 



given the expansion for the velocity (4.1), this simplifies to 

G = 2nr 2 (- - d r ) (yj + ug) (6.2) 

\ r / r=R 1 

The ratio of the effective viscosity of the flow to the kinematic viscosity of the fluid 
is equal to G/G where G is the component of the torque due only to the underlying 
circular-Couette flow, ug. 

Typical nonlinear steady fields in the presence of an imposed axial fields are shown in 
figure @. The main difference between the hydrodynamic and the hydromagnetic cases 
is the axial elongation of the Taylor cells. Not much difference is evident between the 
eigenfunctions and the saturated fields other than that the vortex centres move slightly 
outwards and closer together towards the outflow. 



Figure |2j shows experimental results obtained by Donnelly & Ozima (1962) with mer- 
cury and Perspex cylinders. Also shown are torques for axisymmetric calculations with 
their aspect ratio, r] = 0.95, and results of an amplitude expansion calculated by [Tabeling 



(1981) in the narrow gap limit. 

As the Reynolds number is increased there is good agreement between our numerical 
method, Tabeling's amplitude expansion and Donnelly's experiment, until Donnelly's 
results deviate from both ours and Tabeling's calculations. The points plotted in figure 
^| are time averages as significant fluctuations were observed. Tabeling conjectured that 
this is due to the appearance of wavy modes. With Q = the onset of wavy modes is not 
far above the onset of TVF and in simulations of these modes we find a reduced torque. 

Note that if an axial magnetic field is imposed the onset of wavy modes is signifi- 
cantly inhibited. A detailed investigation of the linear and nonlinear aspects of the wavy 
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Figure 1. Hydromagnetic flow in the presence of an imposed axial field. Q — 180, Re\ = 1.57Je c , 
Re c — 619, rj = 0.95, a — 1.24. The rightmost plot demonstrates that the field lines are dragged 
by the in and outflows. 
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Figure 2. Comparison of torques. Experimental results with rj — 0.95. A, Q — 0; +, Q — 180; 
□ , Q = 652. Solid line, our numerical results. Dashed line, Tabeling's expansion about the 
Reynolds number in the narrow gap limit. 



modes will be presented in the next paper of this series. Here it suffices to remark the 
good agreement between our calculations and the experiment in the weakly nonlinear 
axisymmetric regime. 
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7. Discussion 

In conclusion we have developed a formulation of the governing MHD equations of the 
cylindrical Couette geometry, suitable for timestepping in the nonlinear regime. Results 
agree well with experiments. 

Although the equations do not decouple in the linear part, and we must treat mean- 



flows separately, the formulation is similar to that used by Glatzmaier (1984) in spherical 
geometry. We use potentials for the velocity yet do not eliminate the pressure. This has 
several advantages. Our motivation for adopting such a formulation is that the magnetic 
field then shares the same formulation as the velocity, dramatically reducing the potential 
for error. Only a relatively small part of our code is dedicated entirely to the magnetic 
field; this feature is important for testing, as there are fewer results against which to 
compare our results. Furthermore, it can also accomodate the small Prandtl number 
limit with only minor adjustments. The choice of governing equations which are only 
second order in r makes the method accurate and matrices easily invertible. This feature 
also enables us to take the same radial truncation for all variables, if we desire, simplifying 
implementation a great deal. 

We have opted to use potentials which ensure divergence-free fields. Primativc variable 



for mulations for time integration of the Navier-Stokes equations, such as Marcus (1984) 



and Quartapelle & Verri (1995) in this geometry, do not in general extend naturally to 
the magnetic field. In particular, they are not well suited to the small Prandtl number 
limit, relevant to liquid metals available in the laboratory. 

For the axisymmetric case the expansion by potentials is essentially the same as that 



used by Barenghi (1991) and Jones (1985) and results appear to be very similar in terms 
of accuracy. Although we have one extra equation, the actual form of our equations is 
simpler because we avoided taking the second curl. 

Our method is second order in time and exhibits good temporal stability, We have 
not encountered the difficulties experienced by Hollcrbach (2000b)] and Rudiger fe Feudeij 



(2000) with three-dimensional potential formulations and no-slip boundaries. Using the 
implict Euler method on the linear terms reduces the method to 0(At/Re, At 2 ). With 
finite magnetic Prandtl numbers the magneto-rotational instability (Rudiger & Zhang, 
2001; Willis & Barenghi, submitted) leads to Reynolds numbers which can be surprisingly 
low and the 0(At/Re) error would dominate. 

Results obtained using our method compare well with existing hydrodynamic literature 
with respect to the nonlinear equilibration of Taylor-vortex flow (Barenghi 1991), the 
onset of wavy modes (Jones 1985) and the wavespeed of wavy modes (Marcus 1984). In 
the presence of a magnetic field the results also compare well for the linear stability of 
cir cular-Couette flo w (Roberts 1964), a nd in the nonlinear range th e amplitude expansion 



of Tabeling (1981) and experiments of Donnelly fc Ozima (1962) 



In further work, we will use this method to analyse nonlinear three-dimensional hy- 
dromagnetic Taylor-Couette flow. 

The authors wish to thank Anvar Shukurov and Wolfgang Dobler for stimulating dis- 
cussions and encouragement during this work, and to the referee for helpful clarifications. 



Appendix A. Decay of the magnetic field 

To derive the decay rate of the magnetic field when u = we use a mixture of analytical 
and numerical methods, different from the numerical technique used to solve the MHD 
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m 1 2 3 

k a£ 



(T) 10.634504 2.525434* 6.259403 11.170948 
(V) 9.613411 10.634504 13.640533 18.474045 

1 (l si ) 14.919861 14.760834 16.619549 
(2 nd ) 20.431404 22.010271 25.611050 

2 (1 st ) 45.851458 45.868780 
(2 nd ) 49.822104 51.678920 

Table 5. Decay of the magnetic field, <r£, for r\ = 0.35, a = 3.13. When k — toroidal and 
poloidal modes separate and their slowest decaying modes are given; a is a redundant parameter 
in this case. Otherwise they couple and the first two modes are given. The dominant mode is 
marked *. 



equations. We express the magnetic field in terms of two scalar potentials 

B = VA(Ti)+VAVA(Pz). (Al) 
Each of T, V are expanded and we seek eigensolutions for the magnetic field of the form 

oo 

A(r,6,z)= A km {r)e-' Tt+ ' l{akz+me \ (A 2) 

k,m— — oo 

where a is the decay rate. Substitution into (|2.3| e) yields the Bessel equation 
1 / m 2 \ 

- d r A km (r) + d rr A km (r) + [dr 2 - — )A km {r)=0, a 2 =o-£-a 2 k 2 , (A3) 
which has solution 

A km (r) = A J km J m (ar) + Al m Y m {ar). (A 4) 

Matching at the two boundaries the conditions (3.6-8) for the magnetic field defines the 
problem 

M(a) [T k J m , T km , VL , V km } T = 0, (A 5) 

for the four unknown coefficients. The quantity M is a 4 x 4 matrix and is real if a is real 
(non-oscillatory decay modes). The slowest decaying eigensolution for B is determined 
by the smallest a such that det M(a) = 0. Results are shown in table [j| 
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